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Abstract 

A master equation approach to the numerical solution of option pricing models is 
developed. The basic idea of the approach is to consider the Black-Scholes equation 
as the macroscopic equation of an underlying mesoscopic stochastic option price 
variable. The dynamics of the latter is constructed and formulated in terms of a 
master equation. The numerical efficiency of the approach is demonstrated by means 
of stochastic simulation of the mesoscopic process for both European and American 
options. 



1 Introduction 



As is well-known the seminal Black-Scholes analysis [1] which leads to the fair 
value of an option is based mainly upon the following assumptions. First the 
stock price, i.e. the underlying, follows a geometric Brownian motion. Second 
a hedge position is formed with a portfolio of short underlying and a long 
position of a number of European options. Then an arbitrage argument leads 
to the renowned Black-Scholes partial differential equation determining the 
value of the option. Of course, for simple cases, e.g. constant interest rate and 
volatility, explicit analytical solutions to the equation are known [2]. However, 
for more involved cases one has to rely upon numerical methods. Here we are 
not concerned with deterministic methods, e.g. finite differences, but rather 
with Monte Carlo methods [3]. 
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The idea behind the canonical Monte Carlo method is to exploit the fact that 
the fair value of an option is given by the present value of the expected payoff 
at expiry under the risk neutral measure. Thus the standard Monte Carlo ap- 
proach is based upon the simulation of a geometric Brownian motion for the 
underlying asset until expiry. Then the payoff is computed and discounted up 
to the current time. By averaging over different realizations of this stochastic 
process the current option price is estimated [4]. In spirit this approach is 
similar to the standard Monte Carlo random walk technique for the solution 
of partial differential equations with boundary value problems [5]. 
Monte Carlo methods offer easy to understand numerical algorithms which 
can be easily applied to quite complicated, e.g. path dependent or correlated 
multi-asset, options [6]. In addition they allow for a straightforward inclusion 
of stochastic terms such as the interest rate or volatility. From a numerical 
point of view they are especially suited to problems with many degrees of 
freedom and the algorithm can be easily run in parallel. 
However Monte Carlo simulations of the Black-Scholes equation are usually 
slower than comparable finite difference solutions of the partial differential 
equation. Another disadvantage of the standard Monte Carlo approach is that 
while the application to European options is straightforward, the valuation 
of American options is more involved. Using a generalisation of the canonical 
Monte Carlo method one has to assure that the early exercise condition stating 
that the option price is always above the current payoff is not violated some- 
time in the future. If this happens and the option value is below the payoff, 
the option is exercised and the option value is given by the payoff function. A 
Monte Carlo simulation based on the dynamics of the underlying asset has to 
keep track of all these points in asset time space which makes the algorithm 
ineffective. Nevertheless more advanced Monte Carlo algorithms are available 
[7,8,9,10]. 

Since the Monte Carlo methods described above are based on a model of the 
underlying dynamics of asset values these approaches could be named micro- 
scopic. In contrast now we are going to follow a different strategy, which we 
would like to name mesoscopic, which has already been applied with success to 
simulations of turbulence in fluid dynamics [11,12,13], to the investigation of 
hydrodynamic fluctuations [14,15] as well as to magnetohydrodynamics [16]. 
The same approach has been shown to lead to fast Monte Carlo algorithms 
for the balance equations of nonequilibrium thermodynamics [17], for the sim- 
ulation of chemical reactions [18] and reaction-diffusion processes [19]. 
The basic idea of the mesoscopic approach [20] is to regard the value of the op- 
tion, say V as the expectation value of an underlying fluctuating (mesoscopic) 
stochastic option value, say 6. A master equation for 6 is easily constructed in 
such a way that the expectation value of the multivariate stochastic variable 
9 satisfies the macroscopic Black-Scholes equation. The stochastic process de- 
fined in terms of a master equation is easily simulated and allows for a Monte 
Carlo algorithm for the direct simulation of the Black-Scholes equation. 
This paper is structured as follows. In section 2 the mesoscopic master equation 
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approach is motivated. The dynamics of the underlying asset is modelled as a 
piecewise deterministic process (PDP). The derivation of the fair option price 
for a mesoscopic stochastic option value, which parallels the Black-Scholes 
one, shows that the expectation value of the stochastic option price satisfies 
the "macroscopic" Black-Scholes equation. In section 3 a master equation for 
the mesoscopic fair option value is constructed, and a stochastic simulation 
algorithm is applied to the valuation of European and American options. In 
the last section perspectives and further development of the approach are 
indicated. 



2 Black-Scholes Equation 

The standard textbook derivation of the Black-Scholes equation assumes a 
geometric Brownian motion for the underlying asset which is described by the 
stochastic differential equation 



Here s is the asset value, /is and a 2 s 2 denote drift and variance of the random 
walk. Since dW is the increment of a Wiener process the resulting trajectory 
will be continuous. For more realistic models instantaneous jumps of the un- 
derlying, described by additional stochastic terms based on Poisson processes 
[21,22], may be included in (1). These so called jump diffusion models, are 
used e.g. to incorporate the effect of information about stocks which arrives 
at random times [6]. 

2.1 Black-Scholes equation from a piecewise deterministic process (PDP) 

In contrast to the standard derivation of the Black-Scholes equation based on 
(1), this section demonstrates how the derivation can be paralleled by mod- 
elling the underlying in terms of a piecewise deterministic process [23] . This is 
possible since a diffusion process can be represented as the continuous limit of 
an appropriate jump process. This approach, despite describing the same dy- 
namics as the standard approach, has the advantage of clearly demonstrating 
that it is quite obvious to interpret the option value as a stochastic variable 
and thus to interpret the Black-Scholes equation as the macroscopic expecta- 
tion value of a mesoscopic stochastic option value. In addition the formulation 
in terms of PDPs simplifies the inclusion of additional jump processes since 
one has to deal only with one type of stochastic process. Another important 
point is that for such processes one can easily proceed to a master equation 
formulation, for which powerful numerical algorithms exist, as will be shown 



ds = fisdt + as dW . 
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later on. 

Consider now the stochastic differential equation 

ds = sfidt + 5s dN + - 5s dN~ , (2) 

with (dN^) = s 2 a 2 /(2Ss 2 ) dt. The Poisson increments dN ± are either 0, in 
which case the deterministic evolution takes place, or 1, which describes an 
instantaneous jump of size ±5s. Stochastic process of this general form are 
thus also known as piecewise deterministic processes [24]. 
In the limit 5s — > this stochastic differential equation corresponds to a ge- 
ometric Brownian motion with drift s/j, and variance s 2 o~ 2 [25]. Hence in this 
limit the above equation generates the same dynamics as Eq. (1), but is for- 
mulated as a PDP. 

Using the above stochastic differential Eq. (2) instead of Eq. (1) as starting 
point for the classical Black-Scholes derivation [6] one obtains in the limit 
5s — > for the expectation value of the stochastic option price v 



.dv(s.t), s 2 a 2 ,d 2 v(s,t), .dv(s.t), . , 

where r denotes the riskfree interest rate. If the expectation value V(S, t) = 
(v(s,t)) is identified with the macroscopic option price then one obtains the 
Black-Scholes equation for a European option on a non dividend paying asset. 
Hence the traditional Black-Scholes equation can be interpreted as the expec- 
tation value of a stochastic process. Additional not infinitesimal jump pro- 
cesses can now easily be included, but this of course will change the hedging 
strategy [21,22] and the expectation value of the option price will then follow 
equations similar to the Black-Scholes equation. 



2.2 Macroscopic Black-Scholes equation 



The above derivation assumed a non-dividend paying stock, this assumption 
can be easily omitted to obtain [6,26] 



^ = 4^-,,-,,^^^,,. (4) 

Here a and r again denote the volatility and the interest rate and the continu- 
ous dividend is given by q. In order to numerically solve this partial differential 
equation the relevant variables are made dimensionless and the time direction 
is reversed, i.e. we have, 

2 

S = \n(S/E), V(S,r)=V(S,t)/E, r = y(T-f). (5) 



4 



In this way one obtains 



dV{S,r) _ d 2 V{S,r) 
dr ~ dS 2 



+ (k g ~ 1) 



dV(S,r) 
dS 



k r V(S,r), 



(6) 



where k q = 2 (r — q)/a 2 and k r = 2r/a 2 have been introduced. If the param- 
eters of the Black-Scholes equation (4) are assumed to be constant in time, 
then the analytic solution to the above equation is available, see e.g. [6]. 



3 Master Equation formulation of the Black-Scholes Equation 

As seen in section 2.1 the Black-Scholes equation can be interpreted as a 
deterministic equation which governs the expectation value of a stochastic op- 
tion price. In this section we will describe how the stochastic dynamics of the 
option price can be formulated in terms of a master equation. For the sake of 
simplicity we will consider the Black-Scholes equation for a European call in 
the dimensionless form of Eq. (6). The generalisation to a put is straightfor- 
ward and the application to American options is discussed in section 3.6. 

3.1 General Theory 

A general approach in statistical physics is to derive macroscopic equations 
from a known microscopic dynamics. Here the opposite approach is followed: 
given a partial differential equation for a macroscopic observable we construct 
a mesoscopic stochastic process such that the expectation value of the stochas- 
tic variable is governed by the original partial differential equation. In this 
section the general theory needed for the construction of a stochastic process 
underlying a partial differential equation is presented, for a review see [27]. 
In order to define the underlying stochastic process the state space of the sys- 
tem has to be given. To this end the space variable x is discretized such that 
9\ denotes the state of the system at the discrete position x\, A = 1, . . . , n. 
Hence the state space T of the system is given by T = {9 | 9 e W 1 }. 
The stochastic process is defined by the joint probability distribution P = 
P({9\} ,t) giving the probability of finding the values {9\} at time t. The 
time development of the probability density P is given by a master equation 
of the general form 



The operator A is defined in terms of diffeomorphic maps b acting on the 
state space b : T — > T, and corresponding operators acting on functions 



d_ 

di 



P(9,t) = AP(9). 



(7) 
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F : T — > R in the following way 

B b F(9) = F(b-\9)). 



(8) 



The most general operator A needed in this paper is of the form 

A = cJ2lVet(b ll )]- 1 B b ^-I, (9) 

where Det(6 M ) denotes the determinant of the Jacobian matrix of the map 
and I denotes the identity. In order to generate a valid stochastic process the 
operator A has to fulfil (i) all transition probabilities w are positive and (ii) 
(AF) = in order to preserve the normalisation of the probability density. 
With these definitions the macroscopic equation of motion for a general ob- 
servable F can now easily be computed, one obtains 

^- t (F)^c^(F(b,(e))-F(e)). (io) 

The above equation enables the computation of the macroscopic expectation 
value of a general observable F given an arbitrary multivariate Markov process 
defined by the time evolution operator A. 

For a master equation formulation of the stochastic process one needs the 
transition probability 

w(6, 6) =0^8(6,^(6)). (11) 
Hence the canonical form of the master equation is obtained 

-P(M= J D9~{w(9,9)P(9,t)-w(9,9)P(9,t)}. (12) 

Now the approach sketched above is applied to the Black-Scholes equation. 
Hence the interesting price range of the underlying x = S is divided in n 
discrete points x\, A = 1, . . . , n with distances Aa = x\ — x\-i, A = 2, . . . , n. 

3.2 Stochastic process with uniform discretization 

The easiest possibility to construct the stochastic process is to use a uniform 
discretization of the state space Aa = 51. The Black-Scholes equation (6) 
consists of a diffusive, a convective and a part corresponding to a chemical re- 
action. The stochastic process underlying the Black-Scholes equation will also 
consist of three parts corresponding to these terms in the partial differential 
equation. The total time evolution operator A is thus defined by 

-4bs = Aiiff + A conv ~\~ Achem- (1^) 
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The above time evolution operators are defined through their corresponding 
maps, see Eq. (9). The map 



< = V; Ze'lae- (14) 



describes the diffusive part, while the maps and 

correspond to the convective part and to the chemical reaction term. Hence 
the time evolution operators for the three parts of the stochastic process are 
given by 



A ™=^—aW +A ^- 21 > (16) 
a ol „ 1 — a 

Ahem = ~ J-— C V ~ 1 ( 18 ) 

a ^7 1 — a 

Here the operators A^, B^^C^ are defined according to the general definition 
(8) through their maps a^, b^, c M . In order to prove that the expectation value 
of the stochastic process whose generator is given by (13) really solves the 
Black-Scholes equation one uses Eq. (10) and the projection operator F x (9) = 
9\. One then immediately obtains 



E (0x) ^ + (k - 1) k(9 x ), (19) 

which in the continuum limit 51 — > converges towards the dimensionless 
Black-Scholes equation (6). 

With these definitions the transition probability w{Q, 9) becomes 



{j, n 
+^8(9, c(9)). (20) 

Summarising the Black-Scholes Eq. (6) can be numerically solved by simulat- 
ing the stochastic process 
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Fig. 1. Time evolution of the analytic solution (continuous line) and results of a 
direct stochastic simulation (squares) according to Eq. (13) of the Black-Scholes 
equation for an European call with Exercise price 60, a = 0.2, r = 0.06 for five 
different points in time. The solution of the master equation was estimated by 
averaging 10 realization with a = 0.005 at times T = 10, 7.5, 5, 2.5, 0.05. 

^P(9,t) = A BS P(9) = I D9{w(9,9)P(9,t)-w(9,9)P(9,t)}. (21) 

This is usually done with the help of the following algorithm: 

• initialise 9\ at t = 

• while t < t cn d 

• compute random time step r until the next jump occurs from an expo- 
nential distribution with mean value (r) = 1/iUtot 

■ apply one of the transitions a^,fe M ,c M selected according to their relative 
probability 

By repeating the above algorithm different realizations of 9 are generated and 
thus (9\) can be estimated from a sample of realizations. 
This stochastic process has some interesting features. The parameter a can 
be used to control the size of the fluctuations. In earlier applications of the 
general theory presented above to balance equations of non-equilibrium sta- 
tistical mechanics, the parameter a could be interpreted as the temperature 
of the physical system [27,28]. The larger the parameter a the larger are the 
fluctuations, hence it is intuitively clear that, e.g. in a thermodynamic setting 
a plays essentially the role of the temperature. 

From a numerical point of view the total transition probability w tot is con- 
stant. Thus the error made by taking a constant time step r = l/w tot in the 
numerical simulation vanishes as 0(St). This significantly reduces the number 
of random variables needed. In addition the transition probabilities do not 
depend on the current state of the system. This makes the random selection 
of a transition very efficient, otherwise algorithms as discussed in [29] have to 
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be used. 

Fig. 1 shows a typical result of a simulation averaged over 10 realizations of 
the stochastic process with a = 0.005. Naturally the simulation results fluctu- 
ate around the analytic solution, but these fluctuations are very small if one 
takes into account that only 10 realizations are averaged. The reason for this 
of course is the parameter a. The smaller a, the smaller are the fluctuations 
of the process and the smaller is the number of realizations needed. 
One drawback of this stochastic process is, that it is restricted to a uniform 
discretisation of the dimensionless variable S from Eq. (6). As can be seen in 
Fig. 1 after transforming back to the original variables V and S this leads to an 
exponential distribution of the discrete values of the underlying, see Eqs. (5). 
Hence most of the points are at the left border of the integration region. Since 
the total transition rate performs as 1/51 3 the length of the time step dur- 
ing the simulation and thus the overall computing time depends strongly on 
the discretisation of the underlying. A numerically more efficient algorithm 
thus has to use a uniform discretization of the underlying S which hence re- 
quires much less points. This results in a non-uniform discretization of the 
dimensionless variable S in Eq. (6). 



3.3 Stochastic process with n on- uniform discretization 



To enable a non-uniform discretization of the Black-Scholes equation (6) one 
uses an arbitrary distance of two neighbouring points. The time evolution 
operator ^4 c h em from the previous section does not depend on the discretization 
and is thus not altered. But for the diffusive and convective part of the Black- 
Scholes equation new stochastic processes A^ and A^ v taking into account 
the non-uniform (n.u.) discretisation have to be defined. With the definitions 



°» e » ~ "1 A M _i A 

ft — s. f) -L- m 2fl M A M+ i 

«VH -> tVl-1 + «1 (A M+1 +A M )A M A M+1 



; (22) 



V=P "^" ais ^ 2fjA 2 , (23) 

A 4 I n .a _i_ 2^ M A M _2 ' v I 

for the diffusive and 



" - Vl + «2 A^T 

for the convective part one obtains the two new time evolution operators 
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Fig. 2. Three exemplary realizations of the stochastic process from Eq. (33) (squares) 
and the corresponding exact solution (continuous line) of the Black-Scholes equation 
with Exercise price 60, a = 0.2, r = 0.06 and T = 10. 

' -21, (25) 



Det(r+) " Det(r-) " 



-5 U - 1 



h — 1 

J"-u- __ \ " ; 

The total time evolution operator A^s is now given by 

An.u. _ An.u. _i_ yin.u. , a 
-^BS ~~ -^Miff "T -^conv "^chcm- 



(26) 



(27) 



In order to prove that the expectation value (9\) of the stochastic process 
above really solves the Black-Scholes equation one again introduces the pro- 
jection operator F x and makes use of the general theorem (10) to obtain 



2A A _! 
° In \ A A +A A _! 



) + 



2 A, 



A A +A> 



-2( 



A A _iA; 



+ (*-!)- 



7A+1 



)- 



-k{0 X ). 



(28) 



In the continuum limit the expectation value (0\) of this stochastic process 
thus again solves the Black-Scholes equation, but now on a non-uniform grid. 
The transition probability w{6, 9) becomes 



«W 9) = ^ E W + 6(9, r-{9)) + *-l 2 W «(*)) 



+£e*(M*))- 



«2 



(29) 
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Fig. 3. Time evolution of the analytic solution (line) and results of a direct stochastic 
simulation (squares) according to Eq. (33) averaged over 200 realizations of the 
Black-Scholes equation for an European call with Exercise price 60, a = 0.2, r = 0.06 
at the time points T = 10, 7.5, 5, 2.5, 0.05. 

The parameters ctj are chosen such that the relative size of the transitions is 
smaller than 1, hence they fulfil a < 1, a^/A^ < 1 and 2aq/A^ < 1 for every 
possible transition /i. The smallest ccj, which is aq, determines the scaling of 
the total transition probability 



2n ,, „ . n k 
w = — + (k- 1)— + -. 

aq «2 a 



(30) 



Since aq ~ O(A^) and the total transition probability scales as 0(n/a\) 
the numerical effort to simulate the stochastic process scales according to 
0(n/A^). This is the same dependence on the discretization as in the previous 
section with a stochastic process on a uniform grid. But since now the grid 
can be chosen appropriate one needs much less grid points and the algorithm 
is much faster. 



3.4 Fast stochastic process with non-uniform discretization 



The use of a non-uniform grid reduces the number of grid points to be used in 
a numerical simulation. But there is still a 0(n/ Aj;) 0(1/A^) dependence 
of the total transition rate which makes the algorithm slow. To get rid of this 
dependence the stochastic process corresponding to the diffusive part is again 
improved. In analogy to equilibrium Monte Carlo simulation of lattice systems 
[30,31], where one generates a new configuration of the whole lattice in a single 
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Fig. 4. Numerical root mean square error e of the solution of the Black-Scholes 
equation for an European call with Exercise price 60, a = 0.2, r = 0.06 according 
to the master equation (33) for three different ratios a/N. The parameter values 
for the first simulation (a/N) are ai/N = 10" 7 , a 2 /N = 10~ 5 and a 3 /N = 10" 4 . 



sweep, one defines 



a sweep _ ' 



7"+ j 7"- — 2 1 

Det(t+) Det(t-) 



(31) 



with the map ^ given by 



(32) 



This process thus updates every 9\ at once. The total time evolution operator 
is now given by 



a sweep _ a sweep jn.a. , A 
•^BS ~~ -^diff ' ^conv ' ^chem- 



(33) 



For the transition probability one obtains 



w(9, 9) = — 5(9, t + (9)) + — 5(9, t~(9)) 

+ — E m *$)) + -E W c(9)) . (34) 

Since a± again scales as 0(1/ A^) the total transition probability now also 
scales according to 0(1/A^). This of course makes the algorithm significantly 
faster. 
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3. 5 Analysis of the algorithm 



Fig. 2 shows three exemplary solutions of the stochastic process defined in 
(33) which corresponds to a non-uniform discretisation of the dimensionless 
underlying S. In addition to the expected random fluctuations around the 
exact solution, single realizations also show systematic deviations from the 
exact solution for a wide range of values of the underlying asset. These sys- 
tematic and random fluctuations can be traced back to the two different types 
of stochastic processes entering Eq. (33). The generators -4-conv an d -Achem de- 
scribe single jump processes, where the jump process changes only one option 
value in asset time space, this of course causes the random fluctuations. In 
contrast the systematic deviations stem from -4.^ cp which describes a jump 
process changing all option values at once. 

The numerical estimation of the option price from a sample of realizations of 
the stochastic process is of course not affected by this behaviour of single re- 
alizations. By averaging about several realizations the estimated option price 
converges against the analytic solution, this is shown in Fig. 3. 
Concerning the numerical root mean square error error 



of the Monte Carlo simulation of the Black-Scholes equations one has to dis- 
tinguish two types of parameters entering the numerical algorithm derived in 
the previous section: the number of samples N used to estimate the solution 
V and the parameters ai which control the size of the fluctuations. Figure 4 
shows the root mean squared error e of a Monte Carlo simulation for three 
different values of the ratio ai/N. Different simulations of the same ratio were 
obtained e.g. by decreasing N and all three parameters aj by the same factor. 
This figure clearly shows that the error only depends on this ratio a^/N. In 
addition it can be seen that the error scales with the square root of this ratio. 
This of course is a dependence which is expected for a Monte Carlo simulation. 
In order to assess the numerical performance of the proposed Monte Carlo 
method it is compared with the standard Monte Carlo approach based on 
the simulation of the dynamics of the underlying asset. The numerical error 
in a Monte Carlo solution has two sources. The systematic part of the error 
stems from the finite discretization of the stochastic differential equation in 
the case of the simulation of the geometric Brownian motion and from the 
finite discretization of the state space when one solves the master equation. 
In addition there is a random error from the averaging over different solu- 
tions which decays in both cases as iV _a5 (if a is fixed in the master equation 
formulation). Figure 5 shows the time required to achieve a given precision 
for the two methods in a parameter range where the systematic part of the 
numerical errors are comparable. As can be easily seen the solution of the 




(35) 
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Fig. 5. Comparison of the computation time in seconds needed for the pricing of 
an option by simulation of the stochastic dynamics of the underlying asset and by 
a master equation formulation of the Black-Scholes equation. Shown is the compu- 
tation time vs. the numerical root mean square error e. 

master equation is significantly faster than the standard approach. 
3.6 American options 

Up to now only European options have been considered. How can the concept 
of mesoscopic Monte Carlo simulations be generalised for American options? 
A straightforward generalisation of the algorithm simulating the time evolu- 
tion of the underlying asset is quite involved, since the early exercise condition 
has to be fulfilled. Whenever the option price falls below the payoff function 
the option is exercised and the option value is thus given by the payoff. Hence 
in order to assure that the early exercise condition is not violated sometime 
in the future a straightforward generalisation of this Monte Carlo simulation 
would have to keep track of all these points in asset time space which would 
make the algorithm very ineffective. Of course, there are more advanced tech- 
niques to generalise the Monte Carlo approach to American options, see [10] 
for an overview. 

The proposed mesoscopic approach based on a master equation whose ex- 
pectation value solves the original Black-Scholes equation can be generalised 
to the valuation of American options. Since the time direction in the dimen- 
sionless Eq. (6) has been reversed and the payoff function is used as initial 
condition the generalisation to American options is straightforward. One just 
has to simulate one time step according to the master equation (33) and aver- 
age about the number of samples used to estimate the option price. Wherever 
this option price is below the payoff function it is replaced with the payoff 
and the next time step is simulated. Hence one assures that the early exercise 
condition is fulfilled. 

For this approach to work it is critical that first the samples are averaged 
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Fig. 6. Time evolution for an American call with Exercise price 60, interest rate 
0.07 and dividend yield 0.10 according to the stochastic simulation (squares) and 
the solution of the corresponding partial differential equation (continuous line) for 
the time points T = 10, 7.5, 5. 

and that the early exercise condition is applied thereafter. Applying the early 
exercise condition before averaging would introduce a bias towards higher op- 
tion prices. Every time the stochastic realization is by chance below the payoff 
function the option value is increased. But since the option value is never 
decreased this way it is intuitively clear that one would obtain higher option 
prices. Fig. 6 shows the result of such a simulation for an American call with 
exercise price 60, interest rate 0.07, volatility 0.2 and a continuous dividend 
yield of 0.10. 



4 Conclusions 

In real markets the elegance of the perfect hedge of the Black-Scholes ap- 
proach to option pricing is generally lost. Already a very simple model of the 
dynamics of the underlying stock in terms of piecewise deterministic processes 
shows that the option price itself may be regarded as a stochastic variable. 
However the dynamics of the expectation value of this stochastic option price 
is governed by the Black-Scholes equation. This consideration is the moti- 
vation for the master equation approach to option pricing suggested in this 
paper. 

The essential point of this approach is that the Black-Scholes equation is in- 
terpreted as the macroscopic equation of an underlying mesoscopic stochastic 
process for the stochastic option price variable. By using PDPs one can then 
easily proceed to a master equation formulation of the option pricing prob- 
lem. In contrast to this the usually used microscopic approach describes the 
dynamics of the underlying asset by a stochastic differential equation. 
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The master equation formulation of the option pricing theory offers several 
advantages over the standard approach. This formulation provides a general 
setting in which also other kinds of jump processes may be easily embedded, 
without altering the character of the equations. One may also include addi- 
tional stochastic processes for the volatility and the interest rate and then 
arrive at stochastic volatility and interest rate models [32,33]. This of course 
would lead to different hedging strategies. One possibility to extend the pro- 
posed approach beyond the Black-Scholes equation, which is currently under 
investigation, is to use a hedging strategy proposed in [34,35] 
In addition the master equation formulation also allows for the use of advanced 
simulation algorithms. It was shown that it is possible to construct a master 
equation whose transition probability is constant in time and does not depend 
on the current state of the system. As shown in section 3.1 this allows for fast 
numerical algorithms for the computation of the option price. 
Using the standard Black-Scholes equation as an example we have demon- 
strated how to construct such numerically efficient stochastic processes under- 
lying the partial differential equation. As can be seen from the derivation of 
the general theory in section 3.1, it is clear that this approach is not restricted 
to the standard Black-Scholes equation but can be applied to generalisations 
of the latter. 

The master equation formulation of the Black-Scholes equation is numerically 
about a factor of two faster than the standard Monte Carlo approach. This is 
mainly due to the fact that the proposed mesoscopic approach does not simu- 
late sample trajectories of the underlying asset which, after averaging, results 
in the option price for the initial asset value at a given time. In contrast to this 
the mesoscopic approach works in the whole asset time space and generates 
sample option prices in the whole state space during one realization and is 
hence much faster. 

Another advantage of the master equation formulation of the Black-Scholes 
equation is that it allows for a straightforward generalisation to price Ameri- 
can options in the same framework by just comparing the average option price 
with the payoff function. 

Summarising, the proposed master equation approach is located in between 
the standard Monte Carlo approach for the simulation of the underlying asset 
and the finite difference solution of the partial differential equation. This ap- 
proach tries to combine the numerical efficiency of the solution of the partial 
differential equation with the advantages of the Monte Carlo approach. The 
main advantage of the Monte Carlo approach is that it can easily be applied 
to price options using also advanced stochastic models for the underlying. 
Concluding, the mesoscopic approach we have formulated has been exploited 
to implement a master equation approach to option pricing. Of course, it will 
be of great interest to investigate in further work if features of real markets 
correspond to the additional freedom to choose the size of the fluctuations 
contained in the proposed master equation formulation. This will eventually 
lead to a better understanding of option pricing in real markets. 



16 



References 



[I] F. Black and M. Scholes. The pricing of options and corporate liabilities. 
Journal of Political Economy, 81(3):637-654, 1973. 

[2] J. C. Hull. Options, Futures, and other Derivatives. Prentice Hall International, 
London, 1997. 

[3] Bruno Dupire. Monte Carlo: Methodologies and Applications for Pricing and 
Risk Management. Risk Books, London, 1998. 

[4] P. Boyle. Options: A Monte Carlo Approach. Journal of Financial Economics, 
4:323-338, 1977. 

[5] K. K. Sabelfeld. Monte Carlo Methods in Boundary Value Problems. Springer, 
1991. 

[6] P. Wilmott. Derivatives, The Theory and Practice of Financial Engineering. 
John Wiley & Sons, Chichester, 1998. 

[7] F. A. Longstaff and E. S. Schwartz. Valuing american options by simulation: a 
simple least-squares approach. Rev. Financ. Stud., 14(1):113-147, 2001. 

[8] M. Broadie and P. Glasserman. Pricing American-style Securities using 
Simulation. Technical report, Columbia University, 1997. 

[9] P. Bossaerts. Simulation Estimators of Optimal Early Exercise. Technical 
report, Carnegie Mellon University, 1989. 

[10] P. Boyle, M. Broadie, and P. Glasserman. Monte Carlo Methods for Security 
Pricing. Journal of Economic Dynamics & Control, 21:1267-1321, 1997. 

[II] H.-P. Breuer and F. Petruccione. Burgers's turbulence model as a stochastic 
dynamical system: Master equation and simulation. Phys. Rev. E, 47:1803- 
1814, 1993. 

[12] H.-P. Breuer and F. Petruccione. Stochastic simulations of high-Reynolds- 
number turbulence in two dimensions. Phys. Rev. E, 50:2795-2801, 1994. 

[13] P. Biechele, H.-P. Breuer, and F. Petruccione. Non-equilibrium Monte Carlo 
simulation of a decaying Navier-Stokes turbulence. Phys. Lett. A, 256:147-152, 
1999. 

[14] H.-P. Breuer and F. Petruccione. A Master equation formulation of fluctuating 
hydrodynamics. Physica A, 192:569-588, 1993. 

[15] H.-P. Breuer and F. Petruccione. A master equation description of fluctuating 
hydrodynamics. Physica A, 192:569-588, 1993. 

[16] R. Grauer and C. Marliani. Analytical and numerical approaches to structure 
functions in magnetohydrodynamic turbulence. Physica Scripta, T67:38-42, 
1996. 



17 



[17] H.-P. Breuer, W. Huber, and F. Petruccione. Fast monte carlo algorithm for 
nonequilibrium systems. Phys. Rev. E, 53:4232-4235, 1996. 

[18] H. P. Breuer, J. Honerkamp, and F. Petruccione. A stochastic approach to 
complex chemical reactions. Chem. Phys. Lett., 190:199-201, 1992. 

[19] H. P. Breuer, W. Huber, and F. Petruccione. Fluctuation effects on wave 
propagation in a reaction-diffusion process. Physica D, 73:259-273, 1994. 

[20] N. G. van Kampen. Stochastic Processes in Physics and Chemistry. North- 
Holland, Amsterdam, 1992. 

[21] R. C. Merton. Option Pricing when underlying stock returns are discontinuous. 
Journal of Financial Economics, 3:125-144, 1976. 

[22] J. C. Cox. The valuation of options for alternative stochastic processes. Journal 
of Financial Economics, 3:145-166, 1976. 

[23] M. H. A. Davis. Markov Models and Optimization. Chapman & Hall, London, 
1993. 

[24] H. P. Breuer and F. Petruccione. The Theory of Open Quantum Systems. 
Oxford University Press, Oxford, 2002. 

[25] C. W. Gardiner. Handbook of Stochastic Methods. Springer Verlag, Berlin, 1990. 

[26] P. Wilmott, S. Howison, and J. Dewynne. The Mathematics of Finacial 
Derivaties, A Student Introduction. Cambridge University Press, Cambridge, 
1995. 

[27] H.-P. Breuer and F. Petruccione. How to build master equations for complex 
systems. Continuum Mech. Thermodyn., 7:439-473, 1995. 

[28] H.-P. Breuer and F. Petruccione. Thermostochastics: Heath conduction and 
temperature fluctuations. Physica A, 209:83-95, 1994. 

[29] M. Fricke and J. Schnakenberg. Monte-carlo simulation of an inhomogeneous 
reaction-diffusion system in the biophysics of receptor cells. Z. Phys. B, 83:277- 
284, 1991. 

[30] K. Binder. The Monte Carlo Method in Condensed Matter Physics. Springer 
Verlag, Berlin, 1995. 

[31] D. P. Landau and K. Binder. A Guide to Monte Carlo Simulation in Statistical 
Physics. Cambridge University Press, Cambridge, 2000. 

[32] J. Hull and A. White. The Pricing of Options on Assets with Stochastic 
Volatilities. The Journal of Finance, 42(2):281-300, 1987. 

[33] L. O. Scott. Option Pricing when the Variance Changes Randomly: Theory, 
Estimation, and an Application. Journal of Financial and Quantitative 
Analysis, 22(4):419-438, 1987. 

[34] J. P. Bouchaud and M. Potters. Theory of financial risks. Cambridge University 
Press, Cambridge, 2000. 



18 



[35] J. P. Bouchaud, M. Potters, and D. Sestovic. Hedged Monte-Carlo: low variance 
derivative pricing with objective probabilities. Physica A, 289:517-525, 2001. 



19 



